function derivs = del_log_det_B(hyps, cov_func, n, n_class, X, y, approxF)
% a function evaluating d log.det B d theta_i
% by Mark Norrish, 2011
% hyps: as a vector; [ h_class1 ; ... ; h_class_c ]
% Modifed by ebonilla 09/11/2012
derivs = zeros(size(hyps));

if 0 % comment
  [f, F, bigK, W,  p, P, bigP, sigma_noise, dim, Hyps] = ...
      get_approx_laplace(hyps, cov_func, n, n_class, X, y, approxF);


  dLcond_df = reshape(y, n, n_class) - P; % first derivative of conditional likelihood
  B = bigK*W + eye(n*n_class);


  %disp("B is this asymmetrical");disp(sum(sum(abs(B-B'))));
  % I thought this might have mattered for its partial derivatives
  % apparently not

  BinvK = B \ bigK; % for eqn 5.23
  dL_df = zeros(1,n); % Derivative of the marginal likelihood wrt to fhat

  dK_dtheta = zeros(size(bigK));
  for c = 1 : n_class
    dK_dtheta(:) = 0;   
    crn = 1 + (c-1)*n : c*n; % class-range: i.e. the class-appropriate block
    Bc = B(crn,crn); % consider storing the inverse of this: faster if less stable

      equalities = repmat( (c == (1 : n_class)) , n_class, 1);

      
      % I think these should be computed out of the loop 
      % Mark replies:
      % I don't. W needs be computed forall f_i^c, i.e. n*n_class times, so
      % as is there isn't any extra computation. (Exception: possible vectorisation?)
      % If I took it out, not only would it not save time,
      % but it would waste memory.
      %dW_df = zeros(n_class);
      for i = 1 : n        
	  
	  dW_df2 = 2 * P(i,c) * P(i,:)' * P(i,:) - (equalities + equalities') .* (P(i,:)'*P(i,:)) ...
		      - diag(P(i,:)*P(i,c));

	  dW_df2(c,c) = dW_df2(c,c) + P(i,c);
	  
  %    del_W = zeros(n_class);
  %    for c1 = 1:n_class
  %	for c2 = 1:n_class
  %	  del_W(c1,c2) = -Pi(i,c1)*Pi(i,c)*(c1 == c2) + 2 * Pi(i,c)*Pi(i,c1)*Pi(i,c2) - ((c==c1) + (c==c2)) * Pi(i,c1)*Pi(i,c2);
  %	end
  %    end
	  %for c2 = 1 : n_class
	  %    for c3  = 1 : n_class
	  %       dW_df(c2,c3) =  P(i,c)*(c==c2)*(c==c3) - P(i,c)*P(i,c3)*(c==c2) + ... 
	  %        2*P(i,c)*P(i,c2)*P(i,c3) - P(i,c)*P(i,c2)*(c==c3 + c2==c3);
	  %    end
	  %end

	  % Mark replies:
	  % I'm pretty sure yours is wrong.
	  % Second term, assuming it denotes the first term in my dvi,
	  % should read -P(i,c)*P(i,c3)*(c3==c2)
	  % Last term should read P(i,c2)*P(i,c3)*((c==c2) + (c==c3))
	  % Mine is a transcript of W.m, which gradchek countersigns
	
	  % fprintf('%.6f\n', dW_df(:) - dW_df2(:)); 
	  
	    % disp(sum(sum(abs(dW_df - dW_df')))); % Yeah yours isn't even symmetric.

	  my_range = i : n : n*n_class;
  %disp(abs(sum(sum(BinvK(my_range,my_range) .* dW_df2)) - trace(BinvK(my_range,my_range) * dW_df2))); % gives ~0
	  dL_df(i) = -0.5 * sum(sum(BinvK(my_range,my_range) .* dW_df2));
	  
      end
    
    for hyp_id = 1 : dim
      idx_hyp = (c-1)*dim + hyp_id;
      dKc = feval(cov_func, Hyps(:, c), X, X, hyp_id);
      
      %% Explicit derivative here
      dK_dtheta(crn, crn) = dKc;
      % Mark: question: what? 
      % Why would we want dK/dtheta_h plus dK/dtheta_h2?
      mydkt = zeros(n*n_class); mydkt(crn, crn) = dKc;
      
      %val = (L_B')\(L_B\dK_dtheta);    
      %  val  = solve_chol(L_B',dK_dtheta*W); Surprisingly this still works even if B is non-symmetric
      val = B \ (dK_dtheta*W);
      val2 = B \ (mydkt*W);
      
      explicit =  -0.5 * trace(val2);
      
      df_dtheta = Bc \ dKc * dLcond_df(:, c); % d fhat / d theta_i; by eq 5.24    
      implicit = dL_df * df_dtheta;

      
      derivs(idx_hyp) =  explicit + implicit;
      
      
    end
  end

      

  % B = W*bigK + eye(n*n_class);
  % L  = chol(B)';
  % A = (L'\(L\W));
  % 
  % %% Explicit derivatives here
  % dKc = zeros(n_class*n);
  % for hyp_id = 1 : dim
  %     for c = 1 : n_class
  %         idx   = ((c-1)*n + 1) : c*n;
  %         dKc(idx, idx) = cov_func(Hyps(:, c), X, X, hyp_id);
  %     end
  %     
  %     
  %     for c = 1 : n_class
  %         h = (c-1)*dim+hyp_id;
  %         derivs(h) = -0.5*trace(A*dKc);
  %     end
  % end




  derivs = - derivs; % because it's negloglik not loglik


  return;
else % comment else
  delta = 1e-7;
  %disp(size(hyps))
  for h = 1:prod(size(hyps))
    hyps(h) = hyps(h) + delta;
    lplus = log_det_B(hyps, cov_func, n, n_class, X, y, approxF);
    hyps(h) = hyps(h) - 2*delta;
    lminus = log_det_B(hyps, cov_func, n, n_class, X, y, approxF);
    derivs(h) = (lplus - lminus) / (2 * delta);
    hyps(h) = hyps(h) + delta;
  end
end % end comment

